setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/varscan/bayes_p_errobar")
sel=list.files(pattern = "errobar.txt")
library(ggplot2)
library(dplyr)

for(i in 1:length(sel)){
  file=read.table(sel[i],head=T,sep="\t")
  file1=tidyr::separate(file,normal_rt,into=c("normal_bayes_pvalue","normal_beta0","normal_lower","normal_upper"),sep="_")
  file2=tidyr::separate(file1,tumor_rt,into=c("tumor_bayes_pvalue","tumor_beta0","tumor_lower","tumor_upper"),sep="_")
  file2$total_reads=as.numeric(file2$normal_reads1)+as.numeric(file2$normal_reads2)+as.numeric(file2$tumor_reads1)+as.numeric(file2$tumor_reads2)
  file3=file2[file2$total_reads>=10,]
  file3$normal_bayes_pvalue=as.numeric(file3$normal_bayes_pvalue)
  file3$tumor_bayes_pvalue=as.numeric(file3$tumor_bayes_pvalue)
  file3$unitID=paste(file3$chrom,file3$position,sep=":")
  bias=file3[file3$normal_bayes_pvalue<0.05|file3$tumor_bayes_pvalue<0.05,]
  nobias=file3[file3$normal_bayes_pvalue>0.05&file3$tumor_bayes_pvalue>0.05,]
  bias_normal=data.frame(VAF=bias$normal_var_freq,pvalue=bias$normal_bayes_pvalue,beta=bias$normal_beta0,lower=bias$normal_lower,upper=bias$normal_upper,unitID=bias$unitID)
  bias_tumor=data.frame(VAF=bias$tumor_var_freq,pvalue=bias$tumor_bayes_pvalue,beta=bias$tumor_beta0,lower=bias$tumor_lower,upper=bias$tumor_upper,unitID=bias$unitID)
  bias_normal=bias_normal[bias_normal$pvalue<0.05,]
  bias_tumor=bias_tumor[bias_tumor$pvalue<0.05,]
  biasdata=rbind(bias_normal,bias_tumor)
  biasdata$VAF=as.numeric(gsub("%","",biasdata$VAF))/100
  
  nobias_normal=data.frame(VAF=nobias$normal_var_freq,pvalue=nobias$normal_bayes_pvalue,beta=nobias$normal_beta0,lower=nobias$normal_lower,upper=nobias$normal_upper,unitID=nobias$unitID)
  nobias_tumor=data.frame(VAF=nobias$tumor_var_freq,pvalue=nobias$tumor_bayes_pvalue,beta=nobias$tumor_beta0,lower=nobias$tumor_lower,upper=nobias$tumor_upper,unitID=nobias$unitID)
  nobias_normal=nobias_normal[nobias_normal$pvalue>0.05,]
  nobias_tumor=nobias_tumor[nobias_tumor$pvalue>0.05,]
  nobiasdata=rbind(nobias_normal,nobias_tumor)
  nobiasdata$VAF=as.numeric(gsub("%","",nobiasdata$VAF))/100
  
  datatest=rbind(biasdata,nobiasdata)
  datatest$group="nobias"
  datatest[datatest$pvalue<0.05,]$group="bias"
  p=ggplot(datatest,aes(y=pvalue,x=VAF,color=group))+geom_point(size=1)+scale_color_manual(values = alpha(c('red','#708090')))+theme_light()+geom_hline(yintercept = 0.05,color="red",linetype="dashed")
  fn=paste0(gsub("_bayes_errobar.txt","",sel[i]),"_bias_visualization.pdf")
  ggsave(p,filename = fn,width = 5,height = 4)
}

###可视化的方式2
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/varscan/bayes_p_errobar")
sel=list.files(pattern = "errobar.txt")
library(ggplot2)
i=1
final_bias_tumor=data.frame(matrix(NA,1,6))
colnames(final_bias_tumor)=c("VAF","pvalue","beta","lower","upper","unitID")
final_bias_tumor=final_bias_tumor[-1,]
final_bias_normal=final_bias_tumor
for(i in 1:length(sel)){
  file=read.table(sel[i],head=T,sep="\t")
  file1=tidyr::separate(file,normal_rt,into=c("normal_bayes_pvalue","normal_beta0","normal_lower","normal_upper"),sep="_")
  file2=tidyr::separate(file1,tumor_rt,into=c("tumor_bayes_pvalue","tumor_beta0","tumor_lower","tumor_upper"),sep="_")
  file2$total_reads=as.numeric(file2$normal_reads1)+as.numeric(file2$normal_reads2)+as.numeric(file2$tumor_reads1)+as.numeric(file2$tumor_reads2)
  file3=file2[file2$total_reads>=10,]
  file3$normal_bayes_pvalue=as.numeric(file3$normal_bayes_pvalue)
  file3$tumor_bayes_pvalue=as.numeric(file3$tumor_bayes_pvalue)
  file3$unitID=paste(file3$chrom,file3$position,sep=":")
  bias=file3[file3$normal_bayes_pvalue<0.05|file3$tumor_bayes_pvalue<0.05,]
  bias_normal=data.frame(VAF=bias$normal_var_freq,pvalue=bias$normal_bayes_pvalue,beta=bias$normal_beta0,lower=bias$normal_lower,upper=bias$normal_upper,unitID=bias$unitID)
  bias_tumor=data.frame(VAF=bias$tumor_var_freq,pvalue=bias$tumor_bayes_pvalue,beta=bias$tumor_beta0,lower=bias$tumor_lower,upper=bias$tumor_upper,unitID=bias$unitID)
  final_bias_tumor=rbind(final_bias_tumor,bias_tumor)
  final_bias_normal=rbind(final_bias_tumor,bias_normal)
}
case=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
caseid=unique(case$unitID)

final_bias_normal$VAF=as.numeric(gsub("%","",final_bias_normal$VAF))/100
final_bias_tumor$VAF=as.numeric(gsub("%","",final_bias_tumor$VAF))/100
final_bias_normal$beta=as.numeric(final_bias_normal$beta)
final_bias_normal$group="nosig"
final_bias_normal[final_bias_normal$pvalue<0.05,]$group="sig"
final_bias_tumor$beta=as.numeric(final_bias_tumor$beta)
final_bias_tumor$group="nosig"
final_bias_tumor[final_bias_tumor$pvalue<0.05,]$group="sig"
final_bias_tumor2=final_bias_tumor[!duplicated(final_bias_tumor$unitID),]
final_bias_normal2=final_bias_normal[!duplicated(final_bias_normal$unitID),]

normaldata=data.frame(unitID=final_bias_normal[,6],groupnor=final_bias_normal[,7])
tumordata=data.frame(unitID=final_bias_tumor[,6],grouptumr=final_bias_tumor[,7])
testdata=merge(normaldata,tumordata,by="unitID",all=T)
testdata[testdata=="sig"]=1
testdata[testdata=="nosig"]=2
testdata2=testdata[testdata$unitID %in% caseid,]
testdata2$group=ifelse((testdata2$groupnor==testdata2$grouptumr),"notgood","good")

bn2=final_bias_normal2
bn2$lower=as.numeric(bn2$lower)
bn2$upper=as.numeric(bn2$upper)
bn2=bn2[bn2$unitID %in% caseid,]
bn21=arrange(bn2[bn2$group=="nosig",],beta)
bn22=arrange(bn2[bn2$group=="sig",],beta)
write.table(bn21,"bias_normal_nosig.txt",quote = F,row.names = F,sep="\t")
write.table(bn22,"bias_normal_sig.txt",quote = F,row.names = F,sep="\t")
bn2=read.table("bias_normal_nosig.txt",head=T,sep="\t")
bn2=arrange(bn2,beta,group)
bn2$null=""
bn2$num=1:727
ord=data.frame(unitID=bn2$unitID,num=bn2$num)

bt2=final_bias_tumor2
bt2$lower=as.numeric(bt2$lower)
bt2$upper=as.numeric(bt2$upper)
bt2=bt2[bt2$unitID %in% caseid,]
bt2=merge(bt2,ord,by="unitID")
bt2$null=""
ggplot(final_bias_normal,aes(x=beta,y=VAF,color=group))+geom_point(size=1)+
  scale_color_manual(values = alpha(c('#708090','red')))+theme_bw()

p1=ggplot(bn2,aes(y=beta,x=num))+scale_x_continuous(breaks = bn2$num,labels = bn2$null,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=lower,ymax=upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed",size=2)+scale_color_manual(values = alpha(c('#708090','red'),0.5))+coord_flip()+theme_classic(base_size = 15)+
  theme(axis.text.y = element_blank(),panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")

ggplot(bt2,aes(y=beta,x=num))+scale_x_continuous(breaks = bt2$num,labels = bt2$null)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=lower,ymax=upper,color=group),width=0.5)+
  labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed",size=2)+scale_color_manual(values = alpha(c('#708090','red'),0.5))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()

  
  #######可视化的方式3
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/varscan/bayes_p_errobar")
sel=list.files(pattern = "errobar.txt")
library(ggplot2)
i=1
final_bias_tumor=data.frame(matrix(NA,1,6))
colnames(final_bias_tumor)=c("VAF","pvalue","beta","lower","upper","unitID")
final_bias_tumor=final_bias_tumor[-1,]
final_bias_normal=final_bias_tumor
for(i in 1:length(sel)){
  file=read.table(sel[i],head=T,sep="\t")
  file1=tidyr::separate(file,normal_rt,into=c("normal_bayes_pvalue","normal_beta0","normal_lower","normal_upper"),sep="_")
  file2=tidyr::separate(file1,tumor_rt,into=c("tumor_bayes_pvalue","tumor_beta0","tumor_lower","tumor_upper"),sep="_")
  file2$total_reads=as.numeric(file2$normal_reads1)+as.numeric(file2$normal_reads2)+as.numeric(file2$tumor_reads1)+as.numeric(file2$tumor_reads2)
  file3=file2[file2$total_reads>=10,]
  file3$normal_bayes_pvalue=as.numeric(file3$normal_bayes_pvalue)
  file3$tumor_bayes_pvalue=as.numeric(file3$tumor_bayes_pvalue)
  file3$unitID=paste(file3$chrom,file3$position,sep=":")
  file31=data.frame(file3[,1:4],normal_var_freq=file3$normal_var_freq,tumor_var_freq=file3$tumor_var_freq,normal_bayes_pvalue=file3$normal_bayes_pvalue,tumor_bayes_pvalue=file3$tumor_bayes_pvalue,normal_beta0=as.numeric(file3$normal_beta0),tumor_beta0=as.numeric(file3$tumor_beta0),unitID=file3$unitID)
  file31$normal_var_freq=as.numeric(gsub("%","",file31$normal_var_freq))/100
  file31$tumor_var_freq=as.numeric(gsub("%","",file31$tumor_var_freq))/100
  file31$group=ifelse(file31$normal_bayes_pvalue<0.05|file31$tumor_bayes_pvalue<0.05,"sig","nosig")
}

ggplot(file31,aes(normal_beta0,tumor_beta0))+geom_point(aes(color=group))+scale_color_manual(values = alpha(c('gray','red'),0.9))
